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Abstract 

Pattern formation in uniaxial polymeric liquid crystals is studied for different dynamic closure ap- 
proximations. Using the principles of mesoscopic non-equilibrium thermodynamics in a mean-field 
approach, we derive a Fokker-Planck equation for the single-particle non-homogeneous distribution 
function of particle orientations and the evolution equations for the second and fourth order orien- 
tational tensor parameters. Afterwards, two dynamic closure approximations are discussed, one of 
them considering the relaxation of the fourth order orientational parameter and leading to a novel 
expression for the free-energy like function in terms of the scalar order parameter. Considering the 
evolution equation of the density of the system and values of the interaction parameter for which 
isotropic and nematic phases coexist, our analysis predicts that patterns and traveling waves can 
be produced in lyotropic uniaxial nematics even in the absence of external driving. 

PACS numbers: 
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I. INTRODUCTION 



The study of anisotropic fluids is a very active field of research in soft condensed matter- - — 
in which important results have been achieved very recently.—^ One of the main aspects 
of this study is the obtention of free energies able to describe the transition from orienta- 
tional disordered to ordered states in relation to density or temperature conditions.-^ The 
evaluation of this free-energy is a difficult task- in which mean field theories of anisotropic 
molecules in suspensions or pure systems have proven to be very powerful. The two general 
approaches to the determination of the free-energy function can be categorized as equilib- 
rium and dynamic ones, and lead to different nontrivial structural and dynamic properties 
of nematic liquid crystals.— >2^i^£ 

Equilibrium approaches are focused on the direct calculation of the canonical distribution 
function, the partition function and the free-energy by taking into account the general sym- 
metry properties of the system.— ii£i22i22 For uniaxial systems, this implies that the fourth 
order scalar orientational parameter can be expressed in terms of the second order scalar 
orientational parameter.- Recently, a step forward in the equilibrium approach has been 
given in Ref. [19] where, by using the maximum entropy principle based on Gibbs entropy 
postulate,— a very general method for calculating the free-energy was proposed. In contrast, 
dynamic approaches to the problem are based on a closure of the evolution equations for the 
hierarchy of moments of the non-equilibrium distribution function, typically obtained from a 
Fokker-Planck equation.— An excellent analysis of the physical consistency of several 
dynamic closures proposed in the literature is given in Ref. [18]. 

We consider that these two approaches to the problem are complementary since the 
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equilibrium canonical distribution containing the specific interaction model for the system 
may be used to calculate both, the equilibrium free energy and the corresponding Fokker- 
Planck equation for the non-equilibrium case.-^ 

In general, the results obtained by these two approaches are different even if the same 
interaction model is assumed. To take an example, for the Maier-Saupe interaction model 
in the uniaxial case, which introduces a parameter U measuring the degree of coupling of 
one molecule with the average system surrounding it, the equilibrium approach leads to 
an expression of the free energy in which the dependence on U enters through one term 
(see, for instance, Eq. (20) of Ref. [19]). If entropic effects are adequately taken into 
account, other terms of this free energy appear which are independent of the strength of 
interaction.— In contrast, for this model the dynamic approach leads to expressions of a 
function playing the role of a free energy in which the parameter U multiplies several terms 
that cannot be reduced to that of the equilibrium approach (see, for instance, Eqs. (36) and 
(37) of Ref. [18], among others). In fact, the expression of this function depends on the 
way in which the hierarchy of equations is closed. This situation is usual in non-equilibrium 
statistical mechanics when dealing with closure problems. Different approximations may 
have advantages under different physical conditions. 

We may attribute these differences between free energies derived from equilibrium and 
dynamic approaches to the fact that dynamic closures may incorporate the effect of the 
relaxation of, in principle, all the moments of the distribution function. As a consequence 
of this, it is important to check the performance of the obtained results by following general 
criteria, such as those proposed in Ref. [18]. In addition, because the equilibrium approach is 
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general, it seems appropriate to establish a distinction by referring to a free energy when it is 
derived from equilibrium schemes and to a free- energy like function when the corresponding 
quantity is derived from dynamic calculations.- 1 ^ 

Here, we concentrate our analysis on the dynamic closure approach for uniaxial systems 
and its coupling with the evolution equation of system's mass density. This non-homogenous 
problem is less studied in the literature despite its great interest because it allows to analyze 
the physical conditions in which an appropriate perturbation may produce patterns and/or 
traveling structures in the system. A question of great technological interest. 

The existence of these spatial and temporal structures is analyzed by first using, in 
Section II, the Gibbs entropy postulate^ 1 ^ and the rules of mesoscopic non-equilibrium 
thermodynamics^— in order to derive a Fokker-Planck equation governing the evolution of 
the non-homogeneous distribution of particle orientations. In section III, this equation is 
used in turn to derive the evolution equations of the second and fourth order orientational 
tensor parameters in the homogeneous case. Then, the dynamic closure problems addressed 
in Section IV, where we consider a recent approach^ to determine the forth order scalar pa- 
rameter as a function of the second order one and then obtain the corresponding evolution 
equation for the scalar second order parameter. We also consider an approximation of the 
dynamic equations for both scalar orientational order parameters that allows to propose a 
novel closure relation leading to a novel expression for the free-energy like function. We test 
this closure in comparison with those proposed in Refs. [18] and [19] by analyzing its perfor- 
mance under equilibrium and flow conditions. In Section V we derive the non-homogeneous 
evolution equation for system's density and coupling parameter, and for the scalar second 
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order parameter. Using the mentioned closures, we make a linear stability analysis and 
search for the conditions in which patterns and traveling structures appear. Finally, the 
main conclusions are presented in Section VI. 

II. DERIVATION OF THE ORIENTATIONAL FOKKER-PLANCK EQUATION 

We shall consider a system of N elongated uniaxial molecules in a thermal bath. The 
state of each molecule is specified by the position of its center of mass, r, and by the 
unitary orientational degree of freedom, u, associated with its long axis. The dynamics of 
the system can be described in terms of the single-particle distribution function / (r , u, t) 
that is a conserved quantity and evolves in time by following the continuity equation 

^ = -v l (fv i )-R l (fn l ), (i) 

where, Vi and fli are the conjugate velocities to Xi and Ui, respectively. Here, V« is the i-th 
spatial derivative and Ri is the i-th component of the rotational operator 

du k 

Here, e^k represents the Levi-Civitta antisymmetric symbol and summation over repeated 
indexes will be implicit in through out this paper. 

We will adopt a mean-field approach in which interactions of a single molecule with other 
molecules are represented by the mean-field potential U m { that can be a function of both 
position and orientation vectors, i.e. U m f = U m i(f,u). In addition, we will also consider an 
external non-homogeneous time-dependent force characterized by the potential U c that acts 
on the ensemble of molecules and drives it out from equilibrium. Such a force may represent, 
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for instance, the influence of an electric field.— dLMrM. The spatial and temporal variations 
of the total potential field U = U m f + U e are assumed to evolve in space and time-scales large 
as compared to the time-scales of the relaxation of fluctuations. Accordingly, we assume 
the existence of a local-equilibrium probability distribution, / leq (r,u;t), which is defined in 
terms of the previous potential by the canonical relation 



f lc Hf,u,t) 



z leq (T; t) 



exp 



U (r, u, t) 



(2) 



where 

2 lcq ig 

the corresponding partition function, T is the temperature and fcg is the 
Boltzmann constant. The particular expression for / leq can be obtained from equilibrium 
calculat ions .ILiiL2£ 

The explicit form of Eq. ([I]) in terms of / (f, u, t) can be obtained by calculating the rate 
of entropy production of the system during its evolution in time with the help of the Gibbs 
entropy postulate^ 



S 



-k B J J dfdu f (r, u, t) In 



/ (f, u, t) 
(f, u; t) 



5 lcq , 



(3) 



where iS leq is the entropy of the local-equilibrium state. The entropy production rate o = 
d (S — S leq ) /dt, is therefore obtained by calculating the time derivative of Eq. fl2]), and using 
Eq. (CQ). This procedure yields 



° = T ^ 



fjViifVj+tJiRitf Mi 



+ k B 1 1 drduf^lnf^, 



(4) 



where the non-equilibrium chemical potential, \i = \i (r,u,t), has been defined by 



A/i (r, u, t) = k B T In 



/ (r, u, t) 
f leq (r,u,t)_ 



(5) 
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As we mentioned before, changes in time of the local equilibrium distribution can be 
neglected in the limit of slow external perturbations and, consequently, the second term on 
the right hand side of Eq. (J4]) vanishes. In this case, Eq. (J4]) takes the form 

<y = ~ J J dfdu (pr&in + fQiRit^) , (6) 

where we have performed an integration by parts and assumed an infinite system with 
vanishing probability density at the boundaries in r-space as well as periodic conditions in 
w-space. 

The entropy production rate given by Eq. (jH]) has the form of a sum of products of 
generalized currents, Vi and Qi, with generalized forces, Vj/x and Rifi. According to the rules 
of mesoscopic non-equilibrium thermodynamics, we may follow a linear response scheme in 
which currents are proportional to forces.— 1 ^ Thus, the currents V, and fij are given in terms 
of Vj/i and Rifi by the relations 

Vi = -i^V^/i - MijRjfx, (7) 
Qi = -LijRjH - MijVjfi, (8) 

where Kij, Lij and are Onsager coefficients with Mjj = —My. Replacing Eqs. (GO) and 
(IE|) into Eq. (TjQ), using the expressions for the local equilibrium distribution (T5]) and the 
non-equilibrium chemical potential ([5]), we finally obtain a closed differential equation for /: 
df 



at v * 



Ri 



'A, ( Rjf + ^RjU 



(9) 



Here, we have introduced the diffusion tensors Dij = ksTKy and T>ij = ksTLij. For 
anisotropic molecules both rotational and translational diffusion tensors depend on the ori- 
entational degrees of freedom Wj. The translational diffusion tensor can be expressed in terms 
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of the parallel Du and perpendicular D± (to the symmetry axis) coefficients in the forrn^ 

= D\\UiUj + D± (6ij - UiUj) = D5ij + D a (uiUj - ^ij^J , 

where we have defined D = (D\\ + 2D±) /3 and D a = D\\ — D±. 

At mesoscopic level, the FPE Eq. (Q provides the complete dynamical description of 
an ensemble of uniaxial molecules subjected to mean-field and external forces. It coincides 
with the ones derived in Refs. [1] and [14] and therefore shows their compatibility with the 
second law of thermodynamics. A detailed analysis of the relation of this equation with 
its counterpart derived from kinetic theory arguments in different physical situations can be 
found in Refs. [2] and [21]. In addition it is convenient to mention that the present formalism 
can be generalized to consider the effect of the relaxation of the fluctuating velocities (angular 
and translational) of the molecules in the macroscopic relaxation of the system. This may 
lead to non-trivial diffusion effects, especially in the case of flowing systems.— 1 ^ 

III. HIERARCHY OF MACROSCOPIC DYNAMIC EQUATIONS 

Experimental data characterizing the behavior of a system is more frequently obtained 
through the time course of the components of the tensor order parameter and higher order 
moments of the distribution function /, than from the form and behavior of the distribution 
itself. Thus, it is convenient to derive the evolution equations for these moments. For 
simplicity, this will be done for the case when non-homogeneities can be neglected and 
therefore we shall assume that no external fields are applied and that the description can be 
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carried out in terms of the reduced probability density 

g(u,t) = [ dff{r,u,t), (11) 



which obeys the reduced FPE 

(12) 



9g a 
— — = Rj 
dt 



Rj9 + l^f kjU 



Here, we will assume that each molecule evolves in the presence of a Maier-Saupe mean-field 
potential of the form^ 



3 

Wmf= --UkBTuiUjSy, (13) 



where Sij is the second order orientational tensor whose mathematical definition is given in 
Eq. ( JT3|) below. In addition, U is a parameter that measures the degree of coupling of the 
particle with its surroundings and depends on the molecular structure and interactions.— For 
a concentrated solution of rigid rodlike polymers of length L and diameter b, U is proportional 
to pbL, where p is the number of polymers per unit volume.- 

Following Ref. [1], the moments of g are defined through symmetric traceless tensors of 
rank two, four, etc., corresponding to a multipolar expansion of the orientational degree of 
freedom. In the next, the dyad product of a vector, u, with itself n times, will be denoted 
by U( n ) and 1 will denote the unitary tensor of rank two. General tensors will be denoted by 
blackboard bold characters, A, B, etc. Indices notation will be also used when convenient. 
Finally, [A] s will denote the symmetric part of the tensor A. 

The symmetric traceless tensors of rank two, four and six, are respectively defined by the 
following averages over the distribution g,— 

1 
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w 



< 



"(4)-' [U {2) 1] S + |[11]^, 



(15) 



< 




15 



(16) 



Following an usual scheme in studying the mesoscopic dynamics of nematic liquid crys- 
tals,- we shall assume that rotational diffusion is isotropic, that is, T>^ = T>8^ with T> 
a constant. This assumption is consistent with the preaveraging approximation where a 
constant preaveraged rotational diffusion coefficient is introduced.- 1 ^ 1 ^ 

Thus, using the reduced FPE Eq. (1T2|) . and the previous definitions, the explicit evolution 
equations obtained for § and W are 



Equations (ITT)) and (ITS)) are the first two of an infinite hierarchy of coupled differential 
equations for the moments of the orientational distribution. It is important to notice that, in 
general, the evolution equation for the moment of order n has terms containing the moments 
of order n — 2 and n + 2, thus coupling the whole hierarchy of equations. This is a well 
known result that is explained in whole detail in the classical textbook by Kroger where the 
complete hierarchy is explicitly calculated.- 

An important fact to mention is that the characteristic relaxation time for each one of 
the equations in this hierarchy is proportional to T>~ x . As a consequence of this, it is not 
possible to give a general criterium based on time scales to cut the hierarchy of equations 




(17) 




4 

20£>W + — VU (35 [§§] s + 2§ : § [ll] s - 20 [S • S l] s ) 




(18) 
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and search for consistent dynamic closures. Due to this fact, many dynamic closures to 
the problem have been proposed in the literature to cut the hierarchy. The validity and 
physical consistency of the most representative of them has been analyzed in detail in Ref. 
[18]. However, for comparison purposes, let us recall the classical decoupling approximation 
proposed in Refs. [1, 2, 7], where 

(u (4) ) : § = (u(2))(u (2 )> : S. 

In terms of the tensors § and W, the previous relation can be written in the form 

W:§ = -— §--§■§ + §:§§ + — S:S1. (19) 
15 7 21 v ' 

Indeed, by replacing Eq. ( IT9l into Eq. (117j) . we recover the well known closed expression 
for the order parameter of the Doi-Hess model 



d 

— § = -62? 
dt 



1 - j J § - U ( § ■ § - : § 1 



(20) 



This approximation, as the other ones analyzed in Ref. [18], truncates the infinite hierarchy 
of equations for the moments of the distribution to the lowest possible order. An important 
fact to mention here is that the last three terms at the right hand side of the previous 
equation are proportional to the coupling parameter U. As mentioned in the introduction, 
this is a characteristic result of the dynamic closures in uniaxial systems. The corresponding 
free-energy like function that may be obtained from it by setting the time derivative equal 
to zero, will contain only one term independent of U, which arises from the first term in (1201) 
and is related to the entropic force term of the reduced Fokker-Planck equation (j!2p .- 
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IV. DYNAMIC CLOSURES FOR THE EVOLUTION EQUATIONS OF THE 
SCALAR ORDER PARAMETERS 



We shall consider here a uniaxial nematic liquid crystal. For simplicity, it will be assumed 
to be initially oriented by an external field along the unitary vector n, which is called the 
director. In the following, the dynamics of this phase will be described in terms of the scalar 
order parameters S, W and Z, which are respectively defined in terms of n and the order 
parameters tensors §, W and Z, Eqs. (fT4j) - (|T6l) . by- 

3 

S = -riiSijUj = (P 2 (u ■ ft}), (21) 
35 

W = —nirijWijkirikni = (P 4 (u ■ ft)), (22) 
o 

231 

Z = —ninjn k Z ijk i pq nin p n q = (P 6 (u-ft)), (23) 

where P m is the Lengendre polynomial of order m. Notice that S, W, and Z are zero in the 
isotropic limit and 1 in the completely ordered phase. 

In terms of S, W and Z, the following closure relation can be written for Z in terms of S 
and W, 

Z = ^w { [wsr " n [s ' wir + h [s : Wlir } ' (24) 

where the right hand side contains the symmetric traceless part of the product of WS. This 
explicit relation representing Z is indeed a particular case of the general closure equation for 
moments of arbitrary order derived in Ref.— . 

Using Eq. (I24p the evolution equations for the first two scalar moments can be obtained 
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by projecting Eqs. (|14p and (I15p on n. The result of this procedure is 

^ = -20VW + — VUS 2 + — VUSW- — VUSZ. (26) 
dt 7 77 11 v ; 

Notice that an infinite hyerarchy of coupled equations for the scalar moments is obtained, 
which must be closed at a certain level. In the subsequent sections, two different closure 
relations will be discussed. 



A. Parametric closure 

With the purpose of closing the infinite hyerarchy of coupled equations for the scalar 
order parameters, approximate relationships between 5", W and Z can be obtained based on 
the uniaxial orientational distribution of the Maier-Saupe type^ 

g (u) = G exp [au ■ A ■ u] , (27) 

which is a special case of the Bingham distribution.— Here, a is a parameter, A is a symmetric 
traceless matrix and G is the normalization constant. 

Using Eq. (1271) . S, W and Z can be calculated in terms of a, from which the numerical 
relation between these quantities can be obtained. From this analysis it has been shown in 
Ref.— that W and Z can be very well approximated in terms of S by 

W = S(l-(1- S) v ), (28) 

with ^ = 3/5; and 

Z = S 6 , (29) 
14 



respectively. The first expression can be used in Eqs. (125]) and in order to close the 
hyerarchy of equations upto order S, while the latter allows for closing it at order W. 

Equation ([28]) will be referred hereafter as the parametric closure. It has been shown to 
fullfill the constraints imposed by the nematic symmetry and to be correct in the isotropic 
and totally aligned cases. It also yields the following evolution equation for S 



9S „ m 



j)s-jS* + ^us\i-(i-sy) 



(30) 



Dynamic equations of the form of the previous one, have been extensively used in the 
literature of the field since the pioneer works of Doi and Hess in order to identify a function 
A = A (U, S) through the phenomenological relation^^ 

™ - -L— (31) 

at~ os> (31) 

such that A plays the role of a free energy. As we mentioned previously, it is convenient to 
refer to this quantity as a free- energy like function since it differs from those calculated fol- 
lowing equilibrium techniques. In the present case, the dynamic closure was implemented by 
calculating the equilibrium averages of the second and fourth order orientational parameters, 
and using an interpolation function directly in the evolution equation for the second order 
orientational parameter.— In the following subsection, we will show how the dynamic closure 
approximation is modified when the relaxation of the fourth order orientational parameter 
is also taken into account. 

For the parametric model, A has the explicit form originally obtained in Ref.— 



A= l( 1 _U\ S * + !Ls3+ 12t/ (1 {2 + g (! + »/) [2 + S(2 + u)]} 
2 V 5 J 15 35 (1 + v) (2 + u) (3 + u) 



(32) 
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and predicts that for U < = 4.48, A has only one minimum at S = 0, corresponding to 
the isotropic phase. For U* < U < C/| = 5, A has two minima, one at S = and the other 
one at S > 0, indicating that the system can be found in the nematic or the isotropic phase 
depending on the initial value of S. For U > ?7|, A has a local maximum at S — 0, and the 
nematic phase is the stable one. 

The equilibrium values of the order parameter can be found from the condition 



dA 
dS 



0, (33) 



from which it follows that S eq ({/j* ) = 0.31. 

It should be stressed that although the function A given by Eq. ( 132|) exhibits the typical 
behavior used for describing the I-NPT, it is not completely consistent as long as it is finite 
at S = 1, thus allowing for S to take nonphysical values, i.e. S > 1 for finite energies. 
This problem has been successfully solved in Ref.— , where a thermodynamic free-energy has 
been derived by maximizing the Gibbs entropy postulate ([3]) with the isotropic state as the 
reference one. The obtained expression ensures the constraint S < 1 for arbitrary energies. 

Concerning the high order closure relation ( 1291) . it should be mentioned that when it 
is replaced into Eqs. (1251) and (1261) . and the resulting closed system is used to find the 
equilibrium values of the order parameter S by imposing static conditions, i.e. dS/dt = 
dW/dt = 0, then nonphysical behavior is observed in A since S eq turns out to be a decreasing 
function of U for relative large values of this quantity. Therefore, it should be remarked that 
although Eq. (125]) can be used to approximate very well the values of Z, it can not be 
used consistently in the dynamic equations for S and W by imposing the static conditions 
already mentioned. This suggests that having expressions for higher order moments which 
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approximate very well their exact values, might not be sufficient to fullfill the requirements 
imposed by the dynamic equations. 

B. Dynamic closure 

Here we shall introduce an alternative closure relation for W in terms of S, which has the 
explicit dependence 

W = S%, (34) 

and will be referred hereafter as the dynamic closure relation since it is motivated by the 
explicit form of the evolution equations for S and W, Eqs. ( 125]) and (|26|) . Indeed, it can be 
verified that Eq. (|34|) is exact in two important limiting cases, namely: for U < 1, i.e. in 
the isotropic phase; and for U 3> 1, i.e. close to the completely aligned phase. 

Consequently, this new closure relation is not intended to be exact, but it is an interpo- 
lating expression which is consistent with the time evolution of the scalar order parameters. 
Furthermore, it turns out to approximate the correct parametric relation within the same 
accuracy degree than Eq. ( |28l) . This is explicitely shown in the left panel of Fig. HJ where the 
deviations from the closure relations Eqs. ( 128]) and ( 134]) from the exact numerical parametric 
value of W are presented. It can be noticed that the maximum deviation observed for the 
parametric closure, Eq. (1281) . is 0.052; while the corresponding value for the dynamic approx- 
imation, Eq. (134]) . is 0.065. Notice also that the parametric closure exhibits its maximum 
deviations in the region S ~ 1, where the dynamic closure approximates W better. 
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FIG. 1: Left panel: Performances of the parametric closure, Eq. (|28p (dashed line), and the dynamic 
closure, Eq. (|34p (solid line). AW represents the difference of the approximated expressions (|28p and 
(|34|) , with respect to the numerical value of the order parameter W obtained from the orientational 
distribution (|27p . Right panel: Equilibrium values of the order parameter, S eq , for the nematic 
phase. The black solid line corresponds to the dynamic closure (|36p . which incorporates the effects 
of the relaxation dynamics of S and W. Orange square symbols were taken from Ref. [14] and 
correspond to the numerical solution of the Onsager model. The red dashed-dotted line corresponds 
to the IHK modeU^ that incorporates entropic effects in the Maier-Saupe potential whereas the 
blue dashed line comes from the (KAC)^ parametric closure (|28p . IKO model^ is represented 
by the purple dashed double-dotted line and finally the green dotted line represents the classical 
solution obtained from Doi-Hess model.— 

The dynamic closure relation, Eq. ( 134"1) . yields the following dynamic equation for S, 
dt [\ 5 J 7 35 



(35) 
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Model 


Equilibrium (IHK) 


Parametric (KAC) 


Dynamic (HMS) 


Doi 


IKO 


U{ 


4.59 


4.48 


4.15 


2.67 


6.22 


U* 2 


5.0 


5.0 


5.0 


3.0 


7.34 




0.39 


0.31 


0.41 


0.25 


0.37 



TABLE I: Parameters characterizing the I-NPT in diverse models. U* represents the smallest 
value of U at which the nematic phase can be observed, while is the strength of the mean field 
interaction at which the isotropic phase becomes unstable. 5 eq (U*) is the equilibrium value of the 
order parameter at the I-NPT. This comparison can be extended by considering Ref.—, where an 
extensive study which includes more models can be found. 

and the corresponding function A — A (U, S) for this approximation reads 

^K) s2 4 s3 +t^ (36) 

which predicts the values of the parameters characterizing the I-NPT {/*, U% and S eq (Ui), 
shown in Table [U where these parameters are compared with those obtained from other 
models including the parametric closure. 

The equilibrium values of the scalar order parameter in the nematic phase for the dynamic 
closure can be found from Eqs. (133]) and fl36l) . The solution is shown on the right panel of 
Figured] as function of the reduced interaction strength U/U^- Figure [1] also shows S eq for 
other models including the classical model of Doi and Hess^; the numerical solution of the 
Onsager excluded volume theory as it appears in Ref. [14] ; the model proposed in Ref.— by 
Ilg, Karling and Otinger (IKO), where a generalized mean-field interaction is introduced; 
and the solution obtained from the parametric KAC closure Eq. (|28|) . It can be observed 
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FIG. 2: Oscillations of two elements of the tensor A in terms of time obtained after numerically 
solving Eq. (|38p . In a) we present the diagonal component ^22? whereas in b) the component ^3 
is shown. The values of the parameters used were T> = 0.3, U= 9, Q = 1.5 and the magnitude of 
the shear rate |Vv| = 1.5. 

that our model predicts a wider range of values of U for the coexistence of the isotropic and 
nematic phases than all the other models. Finally, Figure |3] shows the bifurcation diagram 
obtained from the parametric and the dynamic models. 

C. Performance of the parametric and dynamic closures under flow conditions 

In order to rigorously test the performance of the novel closure approximation Eq. 
we will study its behavior under flow conditions. For simplicity, we will restrict our analysis 
to a comparison with respect to the performance of the parametric closure Eq. (12 8|) . which 
in turn has been compared with several closure models as well as with exact results obtained 
from the numerical solution of the FPE.— Our comparison is justified since in the latter case 
it has shown that Eq. (12 8p is valid with high precision for a wide range of conditions. 

We thus consider an elongated molecule immersed in a fluid with a velocity gradient and 
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FIG. 3: Bifurcation diagram of the scalar order parameter S eq as a function of the interaction 
strength U, as obtained from the parametric and the dynamic closures, Eqs. (j28|) (denoted by 
KAC, dashed line) and (fM|) (denoted by HMS, solid line). We also included the equilibrium closure 
reported in Ref.— (denoted by IHL, dash-dotted line). The blue portion of the lines indicates stable 
states whereas the red portion unstable states. The purple portion indicates a stable region for 
KAC and HMS, and unstable for IHK. 

subject to the potential 14. The FPE corresponding to this problem has been calculated 
by considering the hydrodynamic torques which the fluid exerts on the molecule.-^ This 
calculation yields 



dg 1 * 

—— — — LUjRjQ PR; 

dt y 2 



gRi (ujU k K_ 



jk) 







+ Ri 





k B T 3 



(37) 



where w» = f %fcV 'jVk is the vorticity, Ky = \ (ViVj + VjVi) is the symmetric part of the 
velocity gradient, and P = (Q 2 — 1)/(Q 2 + 1) is the shape factor of an ellipsoid of revolution 
with axis ratio Q. 
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Once again, we will consider only a mean field potential of the Maier-Saupe type, Eq. f[T3"j) . 
and calculate the evolution equation for the second order tensor, A = (u( 2 )) = § + |l, in 
the presence of flow directly from Eq. (1371) . This equation involves the fourth order moment 
B = (u(4)), and explicitly reads 

d 

—A + v-VA = A-L-L-A + P(A-K + K-A)- 2PK : B 

—6V (A — -1 ) + 6VU (A ■ A — A : B) , (38) 



3 

where = | (ViVj — VjWj) is the antisymmetric part of the velocity gradient. 

This equation can be closed by using the consistent closure relationship between A and 

m 

B = a[AA] s -2/3[Al] s -2 7 [ll]% (39) 
where the coefficients a, (3 and 7 depend on the scalar order parameters S and W through 

W 

« = 33, (40) 
o a 2W 3 

p= 3 + 2is-r (41) 

7 ^1 + ^-^-^. (42) 
[ 70 45 18 635 V ; 

Thus, we replace Eq. ( 139]) into Eq. ( 1381) and solve the result numerically for the components 
Aij. We consider both, the parametric and the dynamic closure approximations, given by 
Eqs. (125)1 and (|34|) . respectively, in order to compare their performances. For simplicity, we 
assume a homogeneous aligment tensor A, and restrict our analysis to situations similar to 
those reported in the recent literature. Figure [2] summarizes our results. There we present 
the behavior of two different components of A, obtained for diverse values of the parameters 
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P, T>, U and Vu as indicated in the caption. It can be observed that the dynamic closure 
relation, Eq. (151]) . performs quantitatively very well when compared with the parametric 
closure Eq. (128|) . and consequently can be also used to approximate the exact solution with 
high precision at least for the range of values considered here. 

V. PATTERN FORMATION AS A CONSEQUENCE OF THE COUPLING OF 
THE DENSITY AND THE SCALAR ORDER PARAMETER 

The Fokker-Planck equation derived in Section II and the closure approximations dis- 
cussed in the previous section can be used to show that patterns and traveling waves 
may emerge for these systems by following a mechanism different to those discussed pre- 
viously in the literature, where the effect of an external driving was considered.— ~— These 
non-equilibrium structures may occur when the parameter U takes values in the range 
U* < U < U£, that is, when coexistence of isotropic and nematic phases is possible. 

Essentially, the existence of patterns comes from the fact that the degree of coupling U 
may depend on the number density of molecules p. For a lyotropic liquid crystal we may 
write U = (p/p*) U , with U and p* the characteristic energy and density of the system, 
respectively- This dependence indicates that, when increasing the density of the system one 
also increases the interaction energy and the orientational order. This is the basic mechanism 
leading to the INP-T. The formation of stationary and dynamic patterns comes from the 
fact that the resulting equations for S(f,t) and p(f,t), or similarly U(f,t), constitute a set 
of two coupled equations of the react ion- diffusion type.— 

To show this, let us first obtain a dynamic equation for the coupling parameter U(f,t), 
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introduced in the description through the Maier-Saupe potential, Eq. f[T3~j) . We start by 
deriving the dynamic equation for the reduced probability density p defined by 

p(f,t) = J f{f,u,t)du. (43) 

Integrating the Fokker-Planck equation ([9]) over the solid angle yields the following evolution 
equation for the number density 

^ = DV 2 p - l -DV 2 {US 2 ) - l -DV ■ [S 2 VU{p)] , (44) 

where we have used Eq. ( !T3|) and D has been assumed to be constant. Using the relation 
p = (U/Uq) p* with Uq = 5 in the left hand side of Eq. ( |44p and rearranging terms we finally 
obtain 

?LL = d (l - -S 2 ) V 2 U - —DS [VS ■ VU] - —DUV 2 S 2 . (45) 
at \ p* J p* 2p* 

This equation introduces the effective diffusion coefficient D c g = D (\ — -pS 2 ^j . It is inter- 
esting to notice that the correcting factor implies that when orientational order increases in 
the system the diffusion of the particles decreases. However, it is worth stress that molec- 
ular dynamics simulations predicted that the average diffusion coefficient D may increase 
under this conditions.— This competing interplay may lead to an interesting non-trivial 
(non-monotonic) behavior of the effective diffusion coefficient. 

To complete the description, we may use Eq. to derive a general evolution equation 
for the non-homogeneous order parameter tensor Sij (r,t) defined in Eq. (fT4"|) . Multiplying 
Eq. © by UiUj — aud taking the orientational average of the result, an integration by 
parts and some algebra yield the equation 
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FIG. 4: The eigenvalue Ai of the linearized Eqs. (|45p and (|47p as a function of the wave number 
k in cm -1 for different values of the parameters D in cm 2 s _1 (top); P in s _1 and p* in cm~ 3 
(bottom). The eigenvalue becomes negative for p* larger than a certain critical value p* > p c . This 
means that large densities prevent the transition and the appearance of patterns. 



9Sk 
dt 



V 

SVSn - 2- 

J k B T 



u- 



du 



UiUjU k 



dU 

du k 



+DX7 2 Sij + D a V k Vj(^ (uiUj - ^ij) (u k ui - ^8 kl 

+ T^f Vk \ \ iUj ~ ' 



1 



(46) 



where the symbol 8/ du, represents the gradient operator in w-space. It is worth stressing that 
in the proper limiting situations, Eq. (j4~6]) reduces to diverse dynamic equations appearing in 
literature for the tensor order parameter .-^^ Eq. ( 14"6"|) is a particular case of an equation for 
polydomain nematic liquid crystals under shear stresses used to derive microscopic formulae 
for the Frank-Ericksen elastic coefficients.— As indicated previously, Eqs. ( 145|) and (146]) 
constitute a set of two coupled reaction-diffusion type partial differential equations for the 
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FIG. 5: The eigenvalue Ai of the linearized Eqs. (|45p and (|48p as a function of the wave number 
A; in cm~ 1 for different values of the parameters D in cm 2 s _1 (top); D in s _1 and p* in cm~ 3 
(bottom). As in the previous case, the eigenvalue may become negative for p* larger than a certain 
critical value p* > p c . 

scalar fields U (r , t) and S (f, t) . 

With the aim to simplify the description, let us assume that D a = and consider that the 
effective diffusion coefficient of the order parameter, D, is again a scalar constant quantity. 
As a consequence of these assumptions and using Eqs. (132]) and ( |36l) the evolution equation 
for S(f,t) reduces to 



^ = DV 2 S - 6V 
dt 



(47) 



for the five order expansion of the free-energy like function in the parametric KAC closure,— 
and 



dS 
~dt " 

for the dynamic closure. 



DV 2 S - 6V 



(4* 
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FIG. 6: Spatial structures of 5S for different components of a wave vector with magnitude k = 
Gem -1 . The colored regions correspond to different values of S between and 1. The dark purple 
regions correspond to the isotropic state (S = 0) whereas the red regions correspond to a state 
with higher nematic order (S = 0.63). The lighter regions separating the previous ones correspond 
to a nematic order between and 0.63. 

The existence of stationary patterns can be proved by analyzing Lyapunov's stability 
around the nematic equilibrium state Sq shown by the free-energy like functions fl32l) and 
fl36|) . In order to do this, we have to linearize the system of coupled equations (H5]) and 
Wl\ (for the parametric closure), and (jloT) and (H51) (for the dynamic closure) about So by 
assuming S = So + 5S and U = Uq + 5U.— This procedure allows us to calculate the elements 
of the matrix A defining the resulting linear transformation X = A • X, with X_ = (5S, SU) 
and AT the corresponding time derivative. The corresponding elements of these matrices are 
given in the Appendix A. 

The solution of the perturbative system X = A ■ X_ can be proposed in terms of the 
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combination 

X(f,t) = J2^ a e r ^ +Xk \ (49) 

a 

where X a and are the corresponding eigenvalues and eigenvectors of A. 

By performing a numerical study, it can be shown that both eigenvalues are complex. In 
both cases one of the eigenvalues has a negative real part for all the combinations of values 
of the parameters D, p* and Uq, whereas the other one may have a positive real part for 
certain combinations of the parameters. Figures H] and show the corresponding results 
for the parametric and dynamic closures, respectively. These results imply that for the 
adequate combinations of parameters, i.e. physical conditions, the system presents patterns 
that evolve in time in the form of traveling waves.— The projections of these patterns on 
the x — y plane are illustrated in Figure |6] for a wave vector of magnitude k = Qcm~ l and 
three different combinations of its components. Figure [7] illustrates the propagation of the 
patterns in time, for the same wave number. 

Finally, we have also analyzed the case of the exact IHK equilibrium closure given in 
Ref.— by assuming that the spatial and temporal evolution of the scalar order parameter 
may be determined by an equation similar to Eqs. (I45p and 047p . and according to the 
phenomenological approach. The existence of patterns and traveling waves is also possible 
in that case. 

VI. DISCUSSION AND CONCLUSIONS 

In this article, we calculated a Fokker-Planck equation for the single-particle non- 
homogeneous distribution function of particle orientations by assuming an interaction energy 
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FIG. 7: 5S as a propagating wave with k = 6 cm . Each panel shows the configuration of the 
system for t = 0, 4, 8, 12s. Colors mean the same as in Fig. 5. 

of the Maier-Saupe type, which couples the mesoscopic degrees of freedom with the average 
value of the tensor order parameter S^-. The FPE was derived by using the Gibbs entropy 
postulate ([3]), that is also used for equilibrium situations, and calculating the entropy pro- 
duction of the system during relaxation. The obtained FPE is consistent with previous 
descriptions of the dynamics of liquid crystalline phases.- 1 ^^ 

Afterwards, we used the FPE to derive the first two coupled evolution equations for the 
moments of the distribution, Eqs. ( TPT|) and (fTHj) and truncated the corresponding hierarchy 
by adopting two dynamic closure approaches in which the expression of the scalar fourth 
order parameter is expressed in terms of the scalar second order one. The first parametric 
closure, is equivalent to the exact closure reported in Ref.— whereas the second closure, 
called dynamic, represents a novel closure that takes into account the leading term in the 
relaxation dynamics of the second and fourth scalar order parameters, Eqs. fl25|) and (l26j) . 
This novel closure has a good performance when compared with the parametric closure of 
Ref.— and with the exact equilibrium closure of Ref.—, as shown in Figures 1 and 2, and Table 
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1. A comparative discussion on the performance of the parametric and dynamic closures in 
the presence of shear flow was also done. From it may conclude that the dynamic closure 
relation, Eq. ( ]3"%|) . performs quantitatively very well when compared with the parametric 
closure Eq. (128]) . and consequently can be also used to approximate the exact solution with 
high precision at least for the range of values considered here. 

The non-equilibrium part of our analysis was motivated by the fact that Eqs. (1251) and 
fl26|) are coupled in a similar way to Lotka-Volterra equations and, consequently, they allow 
for the possibility of oscillating behaviors.— We search for these behaviors in Section 5, where 
we derived two coupled partial differential equations of the diffusion-reaction type for the 
scalar second order parameter S and the non-homogeneous degree of coupling U. A linear 
stability analysis of these coupled evolution equations showed that patterns and traveling 
waves are indeed possible non-equilibrium solutions for the closure approximations discussed 
in this paper and even for the case of the equilibrium closure of Ref. [19]. 

In summary, we studied the dynamics of uniaxial nematic systems by using the meso- 
scopic non-equilibrium thermodynamics formalism in the context of a mean-field theory and 
within the approach of dynamic closures for the evolution equations for the orientational 
order parameters. The results emerging from this analysis allows us to state that the non- 
equilibrium structures associated to pattern formation and traveling waves are possible for 
these systems without the influence of an external driving. 
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Appendix A: Explicit form of the linearized equations for SS and 6U 



As mentioned in Sec.V, the linearization of the coupled equations for U and 5* about So 
by assuming S = Sq + 5S and U = Uq + SU allows us to calculate the elements of the matrix 
A that enters in the linear transformation X — A • X. The elements of the corresponding 
Jacobian associated to the parametric closure, Eqs. ( )45|) and ( )47|) are 

A ss = DV 2 - ^-V [(24SlU v 2 - 2A^U Q u - 36S$U v + WS U + 7U - 35)] , (Al) 

A sp = -^V [So {6Sy - 6S$v - 12S 2 v + 5S + 7)] , (A2) 



5D 



A pp = ~— [3^oV[/ • V + 3V5 • VC/ + 2S V 2 U + 2U VS ■ V + U S V 2 + U V 2 S ] , (A3) 



5D 



[3S VS • V + 5 2 V 2 + VS ■ VS + S V 2 S ] . (A4) 



In order to search for patterns, we have to propose solutions of the form X_ = X_(t)e lk ' r with k 
the wave vector of the perturbation. Assuming that in Eqs. (1A1I) - (1A4I) we may approximate 
VSq ~ S /C and VS ■ V ~ (Sq/C) ik, we finally obtain the Fourier transformed relations 



A ss = -k 2 D - ^-V [(24^C/ 0l / 2 - 2ASlU u - 36S 2 U v + 10S U + 7U Q - 35)] , 

K = ~T> [So (6S V - 6S 3 u - 12S 2 u + 5S + 7)] , 



A 



fJS 



5D 

App = rSoUo 

p* 

-k 2 D - 5 4s 2 



ik 6 , 2 



(A5) 
(A6) 
(A7) 
(A8) 



where £ = [D/V] 1 is a characteristic length of the system. 

In similar form, the elements of the corresponding Jacobian associated to the dynamic 
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closure, Eqs. ( l4"o"|) and (jlSj) are 



A ss = DV 2 23 

35 



525, 



10/3 



= -—V 



35 



10S"o - 7) Uq + 35 



40/3 



5 125„ ' - 55 - 7 



, (A9) 
, (A10) 



5D 



A pp = -y- [3^oV[/ • V + 3VS • W + 2S V 2 U + 2U VS ■ V + U S V 2 + [f V 2 5 ] , (All) 

5_D 

A ps = DV 2 - — [35 V5 • V + S 2 V 2 + VS • V5 + S V 2 S ] .(A12) 



The corresponding Fourier transformed relations are in this case 



A, 



-k 2 D - -V 
35 



525, 



10/3 



A SP = -^V 



35 



A 



A 



ps 



— T'S'o^O 



105 - 7J Z7 + 35 

S ( 12S 10/3 - 5S* - 7 
zA; 6 



(A13) 
(A14) 
(A15) 
(A16) 
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